#include "header.h"

int main(int argc, char** argv)
{
  work	  wrk;
  params  ctrl;
  aux	  data;
  int	  Emm, K, Jay, l;
  int 	  nref, d = 0, n = 0;
  __float128 *beta_re, *beta_im, *gamma_re, *gamma_im, *rho, *dhi;

  char str[480], line[480], value[480], param[480], value2[480], *dummy;
  sprintf(str, "%s.cfg", argv[1]);

   
  FILE* fh = fopen(str,"r"); ctrl.d = 0;
  if (fh == NULL) err_msg("Cannot open file: input.cfg");
  ctrl.n = 0; ctrl.d = 0;
  while (fgets(line, 480, fh)!=NULL) {
    sscanf(line, "%s\t%s", param, value);
    if (strcmp(param,"#n_poles=") == 0) ctrl.d = atoi(value);
    if (strcmp(param,"#npoints=") == 0) ctrl.n = atoi(value);
    if (strcmp(param,"#resname=") == 0) sprintf(ctrl.resname,"%s", value);
  }
  fclose(fh);
  d = ctrl.d; n = ctrl.n;
  printf("%d\n",d);

  beta_re  = fftwq_malloc(sizeof(__float128)*d);
  beta_im  = fftwq_malloc(sizeof(__float128)*d);
  gamma_re = fftwq_malloc(sizeof(__float128)*d);
  gamma_im = fftwq_malloc(sizeof(__float128)*d);
  rho = fftwq_malloc(sizeof(__float128)*d);
  dhi = fftwq_malloc(sizeof(__float128)*d);

  FILE* fhdat = fopen("FqCqdebug.txt","r");   
  if (fhdat == NULL) err_msg("Cannot open file: FqCqdebug.txt");

  int jj = 0, ind;
  while (fgets(line, 480, fhdat)!=NULL) {
    if (jj == 0) {
      dummy = fgets(line, 160, fhdat);
      dummy = fgets(line, 160, fhdat);
      dummy = fgets(line, 160, fhdat);
      dummy = fgets(line, 160, fhdat);
      dummy = fgets(line, 160, fhdat);
      dummy = fgets(line, 160, fhdat);
      dummy = fgets(line, 160, fhdat);
      dummy = fgets(line, 160, fhdat);
      dummy = fgets(line, 160, fhdat);
      dummy = fgets(line, 160, fhdat);
      dummy = fgets(line, 160, fhdat);
    }
    if (jj < d) {
      sscanf(line, "%s\t%s\t%s", param, value, value2);
      ind = atoi(param)-1;
      beta_im[ind] = -strtoflt128(value, NULL);
      beta_re[ind] = strtoflt128(value2, NULL);
      //printf("%s", line);
      //printf(" %d   %.32Qe   %.32Qe\n", ind+1, beta_im[ind], beta_re[ind]);  
    }
    if (jj == d) {
      dummy = fgets(line, 160, fhdat);
      //printf("\n\n");
      //dummy = fgets(line, 160, fhdat);
    }
    if (jj > d) {
      sscanf(line, "%s\t%s\t%s", param, value, value2); 
      ind = atoi(param)-1;
      gamma_re[ind] = strtoflt128(value, NULL);
      gamma_im[ind] = strtoflt128(value2, NULL);      
      //printf("%s", line);
      //printf(" %d  %.32Qe   %.32Qe\n\n", ind+1, gamma_re[ind], gamma_im[ind]);  
    }
    jj += 1;
  }
  fclose(fhdat);
  dhi[0] = 1.Q - 0.5Q*(beta_im[0] + beta_im[1]);
  //dhi[0] = (beta_im[0] - beta_im[1]);
  for (int k = 1; k < d-1; k++) {
    dhi[k] = 0.5Q*(beta_im[k-1] - beta_im[k+1]);
  }
  __float128 a,b;
  a = sqrtq((dhi[d-3]*dhi[d-3] - dhi[d-2]*dhi[d-2])/(beta_im[d-3] - beta_im[d-2]));
  b = beta_im[d-2] - dhi[d-2]*dhi[d-2]/(a*a);
  
  dhi[d-1] = beta_im[d-2]-beta_im[d-1];
  dhi[d-1] = a*sqrtq(beta_im[d-1] - b);
  
  //dhi[0] = dhi[1] + ((dhi[1] - dhi[2])*(beta_im[0] - beta_im[1]))/(beta_im[1] - beta_im[2]);

  FILE *fhsample = fopen("sample.out","w");
  dummy = fgets(line, 160, fhsample);
  dummy = fgets(line, 160, fhsample);
  dummy = fgets(line, 160, fhsample);

  FILE *fhnew = fopen(ctrl.resname, "r");
  __float128 NV, y0;
  dummy = fgets(line, 160, fhnew);
  dummy = fgets(line, 160, fhnew);
  sscanf(line, "# N = %s\tL2_Residual = %s\ty0 = %s\n\n", param, value, value2);
  //printf("%s\n%s\n%s\n", line, value, value2);
  NV = strtoflt128(value, NULL);
  y0 = strtoflt128(value2, NULL);
  fclose(fhnew);


  FILE *fhres = fopen("poles.dat","w");
  fprintf(fhres, "# 1. chi 2. gamma_re 3. d_chi 4. rho \n");
  fprintf(fhres, "# N = %d\tL2_residual = %Qe\ty0 = %.32Qe\n\n", n, NV, y0);
  for (int k = 0; k < d; k++) {
    rho[k] = -1.Q*gamma_re[k]/dhi[k];
  }
  rho[0] = rho[1] + ((rho[1] - rho[2])*(beta_im[0] - beta_im[1]))/(beta_im[1] - beta_im[2]);
  for (int k = 0; k < d; k++) fprintf(fhres, "%.32Qe\t%.32Qe\t%.32Qe\t%.32Qe\n", beta_im[k], -gamma_re[k], dhi[k], rho[k]);

  fclose(fhres);
}
